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Abstract 

WPHACT (W W and Higgs Physics with PHACT) is a MC program and 
unweighted event generator which computes all Standard Model processes 
with four fermion in the final state at e^e~ colliders. It is based on an 
helicity amplitude method which allows precise and fast evaluations of the 
matrix elements both for massless and massive fermions. Fermion masses 
for h quarks are exactly taken into account. QED initial state and Coulomb 
corrections are evaluated, while QCD final state corrections are included in 
an approximate formulation. Cuts can be easily introduced and distribu- 
tions for any variable at parton level can be implemented. The contributions 
to the processes of neutral Standard Model or Susy Higgs can be included. 
Anomalous couplings effects for the triple coupling can be computed. An 
interface to hadronization is provided and Jetset can be directly called from 
the program. 



^ Work supported in part by Ministero dell' Universita e della Ricerca Scientifica. 
e-mail: ballestrero , accomando@to.infn. it 



Program Summary 



Title of program: WPHACT 

Program obtainable by: anonymous ftp from ftp.to.infn.it in the directory 
pub /ballestrero. 

Computer: DEC VAX, DEC ALPHA AXP, HP/ APOLLO; Installation: INFN, Sezione 
di Torino, via P. Giuria 1, 10125 Torino, Italy 

Operating system: VMS, OVMS, UNIX 

Programming language used: FORTRAN 77 

Memory required to execute with typical data: ~ 500 KByte 

No. of bits in a word: 32 

Subprograms used: VEGAS[|l|], GAMMLN0, RAN20. 

No. of lines in distributed program: ^ 15600 

Keywords: high energy electron-positron collisions, four-fermion final state, VT-pair 
production, Higgs, Z-pair production, LEP2, NLC, QED corrections, electron structure 
functions. Coulomb corrections, anomalous couplings. 

Nature of physical problem 

The forthcoming experiments at the high energy electron-positron collider LEP2 will 
be mainly concerned with WW physics and Higgs search. The production of two W^s 
will allow the direct study of the triple-boson coupling and a precise measurement of the 
mass of the W. The search for the Higgs is of primary importance for understanding 
the problem of mass generation in the Standard Model (SM). Small deviations from 
the SM will be important for discovering possible new physics. Both WW and Higgs 
production will result in a four-fermion final state. It is therefore mandatory to have 
accurate predictions for all physical processes with a four fermion final state in order to 
have full control on signals as well as backgrounds to the processes of interest. The same 
kind of processes will also play a fundamental role in electron-positron accelerators at 
higher energy which will be able to extend Higgs search at higher values of the mass 
and probe triple gauge boson physics and gauge cancellations. 

Method of solution 

Full tree level matrix elements for all processes are computed by means of subroutines 
which make use of the helicity formalism of ref. The number of Feynman 

diagrams in the various channels varies from 3 to 144. The velocity in computing these 
amplitudes that the above mentioned method allows, becomes therefore essential to 
take exactly into account fermion masses and to obtain high precision in a reasonable 
amount of CPU time. 

Different integration variables for the phase space are employed in order to dispose of 
the peak structure of the resonating diagrams for the different processes. The adaptive 
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routine VEGAS is used for the numerical evaluation of the integrals. 

Distributions can be produced to study the behaviour of any variable of interest. For 
simulation purposes, the program can also be used as an event generator that provides 
unweighted events. An interface to standard hadronization packages and specifically to 
Jetset0 is provided. 

Restrictions on the complexity of the problem 

Only for processes with 6's in the final state the masses of the fermions are accounted 
for. Final state radiation is not implemented. Initial state radiation (ISR) is included 
through Structure Functions and no photons pt is computed. QCD corrections are 
introduced only in an approximate way. 

Typical running time 

The running time strongly depends on the process considered and on the precision 
requested. Some examples are reported in Table 1. For the typical final state n~ Ufj, u d 
with ISR the time per call on AlphaStation 600 5/333 is 6. x 10~^ sec. The longest 
time for call is 6. x 10~^ sec. for bbbb. At Lep2 energies, 5 M calls (about 5 minutes) 
are used to obtain for fi~ u d a. cross section (with ISR) with a typical estimated 
relative error of 2 x 10~^ sec. The same process can be evaluated in about 40 sec. with 
less than 1 M calls at permill level. All above running times have to be multiplied 
approximately by a factor 3 for an AlphaServer 2100 4/200 computer. 

Unusual features of the program: 

REAL*8 and COMPLEX* 16 variables and STRUCTURE declarations are used. Compilation 
on Hp/Apollo stations has to be performed with -k option. 
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Long Write-Up 



1 Introduction 

After LEP 1 precision measurements and top discovery at Tevatron, which have shown 
an impressive agreement of the Standard Model (SM) with data, a new era is starting in 
which the mechanism of spontaneous symmetry breaking and the non abehan structure 
of the model will be directly tested by the experiments. At LEP 2 it will be possible 
to measure the contribution of trilinear gauge boson couplings, and for the first time 
the production of two W's and two Z's will be seen. The properties of the W's will 
be measured with great accuracy and these will also contribute to put more stringent 
limits on the Higgs mass. At the same time, direct searches for the Higgs boson will 
allow to find it, if its mass is not greater than ^ 100 GeV. These studies and the 
search for the Higgs will be extended to LHC and to the next e^e~ collider which will 
presumably approach the TeV range. Theoretical arguments claim for the onset of a 
'new physics' regime at this scale, so that LEP 2 and the future machines will search 
not only for new particles, like the MSSM Higgs, but for any deviation from the SM 
predictions. 

It is evident in this scenario that theoretical predictions must reach a high accuracy 
to confront and analyze the data. WW and Higgs physics and their radiative corrections 
have been thoroughly studied in the past (for a complete review on this subject see 
ref.@). It has however been realized that on shell predictions may not be accurate 
enough, as the measured final states will not correspond to, say, two W's or a Higgs 
and a Z, but rather to 4 fermions. This implies that one has to deal with irreducible 
backgrounds, i.e. contributions to the matrix element for the envisaged final state 
which do not correspond for instance to two W^s production and decay and cannot 
be separated from it, even if they can be reduced with appropriate cuts. Moreover 
many different final states are not experimentally distinguishable, and one has to take 
into account that the properties of the W and Higgs bosons can be only reconstructed 
by a careful analysis of all four fermion final states. For top physics, trilinear and 
quadrilinear coupling studies and for the Higgs, if it will be as heavy as to decay into 
two ly's, also six fermion final states will have to be analyzed at future e+e~ colliders. 

Several codes for four fermion physics have been produced in the last years 
0-[|T^, and they have been used for phenomenological studies mainly concerning WW 
and Higgs Physics [|lT|-|]l2|-|jl3|. Only some of them can produce accurate results for 



all four fermion processes, with the inclusion of all the relative Feynman diagrams. 
During last year's LEP 2 Workshop the results produced by the different codes have 
been extensively compared The codes have very different characteristics. 

Some of them are classified (see reLQ) as semi- analytical [|12| or as deterministic]^]. 



All the others belong to the broad Monte Carlo's (MC) class, where a further subdi- 
vision can be made among unweighted event generators and Monte Carlo integration 
programs. The semi-analytical codes perform as much as possible of the integrations 
analytically, leaving only low dimensional integration to be performed numerically. The 
MC and deterministic programs on the contrary perform the whole of the 9-dimensional 
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integration (including ISR) numerically. This implies that by their own nature the semi- 
analytical codes may easily reach an extreme precision, but they cannot implement all 
cuts: normally only those on the invariant masses are viable. The MC and Determin- 
istic programs can implement all cuts, but of course they are normally slower and less 
precise. The deterministic program[Q implements all cuts analitically, while only few 
of them can be implemented in MC's as limits on the integration variables. 

WPHACT belongs to the family of the MC integrators and event generators. It has 
been developed only in the last year but it has been compared and tested continuously 
both in LEP 2 Workshop 0-[|lOl-[|l4| and in phenomenological studies[^3|]. One of the 
main characteristic of the code is that of using a new helicity amplitude formalism^]- 
0] which allows to compute matrix elements in a very fast and precise way. As a 
consequence, the code can reach high precision in a relatively short time. This is 
particularly useful for example when one wants to use WPHACT to produce distributions 
at parton level. If one requests that the range of a variable be divided in a large 
number of bins, of the order of 100 say, one can reach very low errors on each bin only 
producing a number of weighted events of the order of ten millions. And this is precisely 
what WPHACT can do in a reasonable time (see Table 1) even for the most complicated 
processes as e^e~ —>■ bbbb with all diagrams and taking b masses into account. 

As already mentioned in the program summary, WPHACT allows to compute QED 
initial state and Coulomb corrections, and also QCD final state corrections in a 'naive' 
formulation, exact in the limit in which only double resonant diagrams are considered 
and no cuts are imposed. Possible cuts are explicitly provided and any distribution at 
parton level can be implemented. WPHACT is also an unweighted event generator. The 
events can be produced while evaluating the cross section and also after this has been 
computed. In this second case any number of events can a priori be requested. An 
interface to hadronization is in the code and it may be linked directly to Jetset. As far 
as 'new physics' is concerned. Anomalous Couplings effects for the triple coupling and 
cross sections for neutral MSSM Higgs h and A can be computed. 

In presence of unstable gauge bosons the imaginary part of their propagators violates 
gauge invariance in tree level computations of processes like those we are considering. 
The way to restore it and get a reliable result is that of including relevant parts of 
fermionic corrections, thus fulfilling Ward identities [|1^]. From a numerical point of 
view these corrections become relevant only for some particular processes and cuts. 
For instance in the process e~ u d no discrepancy between approximate and correct 
computation is appreciable if a cut of the order or greater than ^ 5° is applied to 



the angle between the electron and the beam||13[. The present version of WPHACT does 



not include these corrections. Also final state radiation and effects due to transverse 
momenta of QED radiation are not computed at present by WPHACT. Subtle theoretical 
problems are connected to these issues and probably only a full computation 0{a) can 
assess them. As far as final radiation is concerned, this can in any case be introduced 
via Jetset when WPHACT is used as an event generator. 

In the following we will describe the main features of the code and the way to make 
use of it. Some useful examples of test runs will be given at the end. 
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2 General Features 



2.1 Processes 

WPHACT computes all SM processes with four fermion in the final state at e^e~ colliders. 
Final states with t quarks are not considered, as the t's are known to decay immediately 
to 3 other fermions. 

The processes are enumerated in Table 2 and 3. They are divided in 4 classes. The 
first (CC) and the third (NC) contain all processes which have only diagrams with 
charged or neutral currents respectively. In the second (Mix) they have both kind of 
diagrams. In the last class (NC+Higgs) we have included the processes with 6's in the 
final state, where we may have important contributions from the Higgs. These are the 
only processes in which WPHACT accounts for Higgs diagrams and it treats the h quarks 
as massive both in the phase space and in the matrix elements. 

In WPHACT the momenta of the particles in the final state, as reported in the tables, 
are respectively p3, p4, p5, p6. The order of the final particles is important when 
one wants to choose a certain set of cuts for the final state, or to ask the program to 
compute some distribution at parton level. It has to be noticed that for CC and Mixed 
processes the order in which the momenta are passed to the high energy standard 
COMMON/HEP EVT[|l6l and hence also to Jetset may be different from the one used by 
the program: in the common particles 3 and 4 are respectively the particle and the 
antiparticle that correspond to a , 5 and 6 to a W~ . 

Many different final states give the same cross section at parton level. This is the 
case for instance of fi~ d d and t~ s s if the mass differences are neglected. For 
this reason to each final state two fiags are attributed, iproc and ich, which serve to 
identify it. The first refers to the group which has the same diagrams (whose number 
is indicated in the first column) and cross section, the second to the specific final state 
of the group. It has to be noticed that for CC processes charge conjugate final states 
belong to the same iproc. The amplitudes for a CC final state and its charge conjugate 
for a given set of four momenta, are not equal. To obtain the same amplitude one must 
consider the parity transformed final state of the charge conjugate. This implies that 
if the cuts are not invariant under parity transformation (i.e. both a set of final four 
momenta and its parity transformed are accepted or not accepted ) the two cross 
sections will be slightly different. An example of this is given by e~^e~ ^v^ud and 
e^e' fi^Ufj^ud when in both processes the cut angle of the /i with is different 
from that with e". WPHACT accounts exactly also for these slight differences. The main 
reason to specify the exact final state instead of just indicating the group is however 
to distinguish among final states when hadronization via Jetset is performed, as this 
depends on the final state particles. 

The production of two W^s which decay in two fermions each is described by 3 
Feynman diagrams. The cross section corresponding to these 3 diagrams only is often 
referred to as CC3p. Of course CC3 does not correspond to physical processes, and 
the three diagrams alone do not constitute a gauge invariant set, but for some energies 
and final states and cuts they might be a reasonable approximation. For all CC and 
CC+NC final states, WPHACT can compute the complete process or the CC3 contribution 
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to it. 

All neutral current processes are normally computed at order a^. When there are 
four quarks in the final state, there are however contributions O(a^a^) of diagrams in 
which a gluon is exchanged between the two quark lines. This contribution is of course 
enhanced by the coupling and can be relevant for some processes, energies and cuts. 
For instance we have found[0 that for a{e^e^ bbbb) at Ecm = 175 GeV and with 
f^bb ^ 20 GeV, it is more than twice the pure electroweak contribution, while it becomes 
only about one third at Ecm = 192 GeV. WPHACT may or not include this contribution. 
If one includes it, one must be aware that part of it may also be accounted for by 
two quarks final state computations if parton shower and hadronization are added. 
The interplay among these two ways of treating such a contribution is surely strongly 
dependent on cuts used and deserves a careful analysis case by case. 



2.2 Helicity amplitudes 

All amplitudes for the processes computed by WPHACT are evaluated with the helicity 
formalism of ref . , which is based on the insertion in spinor lines of a complete set of 
states for every fermion propagator. These states are eigenstates of ^, where p is the 
momentum fiowing in the fermion propagator. They are chosen to be generalizations of 



the spinors used in ref. . In this way one needs not to decompose every |^ in terms of 
the external momenta, as it is the case for some other helicity method, and the numer- 
ator of the fermion propagator assumes a very simple expression. The computation is 
in this way reduced to the evaluation of one r matrix for every vertex and to combine 
them together. For the definition of these 4x4 matrices we refer again to ref. 0. In 
this base the matrix corresponding to the numerator of every fermion propagator is di- 
agonalized for massless lines and for massive lines the non zero off diagonal elements are 
simply given by the mass itself. It has to be noticed that a r matrix fully describes the 
vertex both in the case of insertion of an external particle and in that of the insertion of 
a whole piece of diagram. Moreover, combining together two r matrices corresponding 
to two vertices, one obtains a new r matrix corresponding to the two vertices and so 
on. These simple facts allow to achieve a great modularity in the computation, to keep 
track of partial results and to set up a recursive scheme which computes and stores for 
later use sub diagrams of increasing size and complexity. Moreover the massive case is 
not more complicated than the massless one. Only more helicity indices are of course 
needed. As a consequence, the codes for massive amplitudes written in this way are 
not much slower, as it is normally the case, than those with massless fermions. 

The code for WPHACT amplitudes has been completely written with the help of 
PHACT (Program for Helicity Amplitudes Calculations with Tau matrices). This 
program implements the method described above in a fast and efficient way. It essen- 
tially writes directly the optimized fortran code for every r matrix corresponding to 
a vertex and for combining r matrices together. It also writes the code for external 
gauge boson polarization vectors, triple couplings, and so on. With it, one avoids as 
much as possible computing expressions which will turn out to be zero. Moreover in 
the computations of the amplitudes we have avoided any call to external subroutines 
and functions which might be time consuming for such a repetitive part of the program. 
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This of course leads to a somewhat longer program than usual. 



2.3 Phase space and integration variables 

To describe the phase space, the four final fermions are divided in two couples. Every 
couple corresponds to two particles that can decay from a in CC, from a Z or a Higgs 
in NC contributions. Take for instance the case in which /4 and /s /g are grouped 



together. Natural variables are then the two invariant masses 77134 = J {p3 + p4) 



'^se = Y (p5 + p6y and the decay angles of one particle for each couple ^3, 03, 6t,, 05 in 
the rest frame of the couple itself. 

Whenever the energy is such that a couple of particles can have an invariant mass m 
equal to M, the mass of a real W or Z or Higgs, the corresponding amplitude squared 
will be almost proportional to a Breit-Wigner peak: 

(1) 



It is therefore convenient to use instead of the invariant mass m of the couple an 
integration variable proportional to : 

- 

X = arctg . (2) 

7 

Every variable is always transformed to the interval -h 1. When substitution is 
performed on both invariant masses, we refer to it as double resonant mapping. If 
only one or no substitution is performed we have respectively a single resonant or non 
resonant mapping. In WPHACT there is the possibility to choose among these mappings. 
Most of the results at LEP 2 energies are better obtained with the double resonant 
phase space. For some processes one might have a peak also for low invariant masses, 
due to photon propagator contributions. Whether or not this peak becomes relevant 
depends strongly on the cuts. When such is the case, one has a variable m to which 
there correspond two peaks: the Breit-Wigner and the photonic one. After having 
tried different mappings, we have found that, given the fact that VEGAS is an adaptive 
routine, it is better just to use a non resonant mapping for m in this case. 

When one has to deal with a Mixed (or NC+Higgs) process, the peaking structure 
of CC (Higgs) and that of NC contributions is different. In such cases WPHACT auto- 
matically integrates separately the two contributions. To one of them it is added the 
interference. 

The contribution to e'*"e~ — > bbbb of diagrams with h A intermediate states and that 
oi h Z give two {A and Z) Breit Wigner peaks on the same variable m. In this case a 
double mapping of the type (|^) is performed simultaneously on it. 

With ISR two more integrations are to be performed and we use also for this case 
a change of variables to level the form of the distribution functions, as we will see in 
the following section. 

All integrals are computed with VEGAS. For this routine it is convenient to use 
normally more than one iteration. In every iteration the integral is evaluated and at 
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the end the various results are combined together. This allows to optimize the number 
of points computed in the various regions of the integration variables. It may also 
be useful to perform some thermalizing iteration with a lower number of points to be 
evaluated. The thermalizing iterations are just used to adjust the grid and not for the 
final result. 

The adaptivity of VEGAS is such that, even if most cuts are implemented in the 
program with the use of if statements which act as 9 functions, usually this does not 
correspond to a sensible lost in time and precision. The cuts on the invariant masses 
which are function of integration variables are implemented directly on integration 
limits. 

2.4 QED and QCD corrections 

QED initial state radiation is taken into account in WPHACT via the structure function 
approach [[l^] . In it, the cross section a{s) is obtained by means of a convolution with 
functions D{x, s) which account for the radiative emission by the initial particles of a 
fraction x of their energy: 

a{s) = jdx.dx,D{x.^s)D{x,^s)a{x.x,s) (3) 

The expression for the 0{q^) D[x, s) in leading log approximation used by WPHACT is: 



D(x, s) 



exp 



2r (1 + 1/5) 



Pil-xf/'-'-^il + x) 



(4) 



+ 



32'^ 



4(1 + x) ln(l — x) H — In X + 5 + X 



with (3 = lajn (ln(s/mg) — 1) and 7^; the Euler constant. 

Given this form of structure functions, it is convenient for the MC to perform a 
change of variables in order to flatten out the peak of the distributions near x = 1. 
Thus we have chosen instead of X\ and X2 the integration variables 



:i - 



1,2) 



(5) 



Another important QED correction to be accounted for is the so called Coulomb sin- 
gularity, which is due to the electromagnetic interaction of the two W bosons at low 
velocity. It regards of course only CC and Mixed processes near the threshold energy 
for production of two W^'s. In that region however its correction amounts to a few 
percent. In WPHACT this correction is computed with the method of ref . [p!9[] . 

Final state QCD corrections are taken into account by means of the so called 'naive' 
approach. In it, the various contributions to the amplitude are multiplied by the 
corrections relative to the decay width of the VT's, of the Z's or of the Higgs. This 
approach disregards all corrections to non resonant diagrams and corrects exactly the 
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decay vertices only if no cuts are applied. Nevertheless they represent at present an 
estimate of a sizeable effect. 

For NC, CC and Mixed processes the QCD corrections appear in an overall factor: 
(1 + ^) for semi-leptonic case and (1 + 2^) for non-leptonic one. Gluon exchange 
diagrams and CC{NC) ® QCD interference are not multiphed by a correction factor 
as we only consider the first order in ctg. For the processes involving h quarks, Higgs 
exchange diagrams receive the following strong correction: (1+5.67^) for semi-leptonic 
case and (1 + 6.67^) for non-leptonic one. The NC ® Higgs interference is multiplied 
by (1-1-3.335^) for semi-leptonic case and (1-1-4.335^) for non-leptonic one. By default 
as in the preceding formulas is taken at Mw scale for CC diagrams, at Mz scale for 
NC and NC-|-Higgs ones. Its value can however be changed in a DATA statement. For 
consistency, when corrections are applied to the vertices they should also be present in 
the same way in the widths appearing in the propagators of the bosons. This is surely 
achieved if one chooses to let WPHACT compute the corresponding total widths. 



2.5 Susy and Anomalous couplings 

Besides SM Higgs processes in two 6's and two other fermions, WPHACT computes also 
Susy neutral Higgs production in the same final channels. The only MSSM process 
which cannot be deduced just changing couphng constants with respect to the SM is 
e+e~ — > hA — > hhhh, where h is the hghtest CP-even and A the CP-odd Susy Higgs. 
The cross section for the above diagram is computed as well as all other contributions 
to 4 fe's final state. As for all other Higgs processes, the contributions of all diagrams 
containing the Higgses can be optionally separated from the rest. The Susy parameters 
to be given in input are the mass of the pseudoscalar Higgs A, and the ratio between the 
two vacuum expectation values tg{l3) . In the simplest version of the MSSM, all Higgs 
masses arc predicted in terms of these two parameters. At one-loop these predictions 
are substantially modified and an additional dependence on the top mass rrit and on 
the common squark mass mi is introduced. We have used the following relations: 

Ml = ^[Ml + M| + e/sin2/3] 

-\ {[{Ml - Ml) cos 2/3 + e/sin2/3]2 + (Mj + M|) W2/3}'^' , (6) 



where 



STT^M^sm^^ \ mj) 

The squark mass scale mi has been chosen to be 1 TeV. The mixing angle a in the 
CP-even sector, which together with (3 determines all couphngs of the MSSM Higgses, 
is defined by 

(M^ + M|)sin2/'j 
(Ml - M|)cos2/3 + e/sm^/? ^ ' 

As far as Anomalous Couplings computations are concerned, we have imple- 
mented those relative to the trilinear vertex. Starting from the most general effective 
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lagrangian|pO||, one gets for them nine possible couplings just imposing Lorentz and 
electromagnetic invariance. These are further reduced to 6 if CP invariance is im- 
posed, and 5 of them are separately C and P invariant. Among the various possible 
parametrizations of these 6 quantities, we have implemented that of ref pi | . 



2.6 Distributions and unweighted event generation 

One of the main purposes of a dedicated code as WPHACT is to perform phenomeno- 
logical studies and confront theoretical predictions with experiments. To this aim, the 
possibility of computing differential cross sections or distributions is extremely impor- 
tant. Special care has been devoted to this aspect and practically any distribution at 
parton level can be computed while the total cross section is evaluated. Explanations 
of how to request and implement distributions are given in section 3.1. In practice, 
one has just to write down in an include file the definition of the variables to be dis- 
tributed in terms of the 4-momenta of the outgoing particles and to specify in input 
the number of bins and the interval for every such variable. Files with the extension 
.dat will contain in output the cross section and the evaluated statistical error relative 
to every single bin. A large amount of bins (of the order of 100) per variable with a 
low statistical error can easily be achieved. 

If the distribution refers to leptons in the final state, it will correspond to some 
directly measurable variable. For quarks, the hadronization process might prevent this 
possibility. Distributions at parton level are however much faster to obtain than the 
ones with particles which have to be computed after parton shower and hadronization. 
They constitute a very effective tool to study the physical problem at hand and to 
deconvolute perturbative effects from non perturbative ones. Examples of two distri- 
butions obtained with WPHACT are given in figs. 1,2. 

With WPHACT it is also possible to generate unweighted events and to store their 
4-momenta. From them one can of course successively produce any distribution. To 
reach low statistical errors in this way a high number of events has to be stored, but 
they can be produced in a reasonable time: for instance after having evaluated the cross 
section, it takes 48 min. to produce 500000 unweighted events for MIX19 e^e~ — > /i^ 

WPHACT can generate unweighted events with the hit or miss method while evalu- 
ating the integral. In that case the first effective iteration after thermalization is used 
to find the maximum and the second to generate unweighted events. One may also 
produce them after having evaluated the integral. This feature can be used both to 
readjust the maximum if some event exceeded it with the first method, or to ask for a 
predetermined number of events to generate. When unweighted events are produced, 
they may be passed directly to Jetset using the routine AB_LU4FRM. So one can for in- 
stance run the first time just to find an appropriate maximum and a correct cross section 
(and eventually to produce parton distributions). Afterwards one verifies whether it 
is convenient to multiply the maximum for an appropriate factor and with a second 
run one generates the desired number of events. To save time, the parton shower and 
hadronization procedure with the link to JETSET may be activated only in this second 
run. 
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As explained in sect. , the integration of Mixed and NC+Higgs complete pro- 
cesses is performed in two steps, with different mappings of the phase space. Also 
the event generation will be performed in two steps in these cases. This implies that 
the events that come first will have been selected with the CC (Higgs) contribution 
and the remaining with the NC contribution. The resulting sample will have the right 
proportion of events, but from them one cannot take away an arbitrary part. If for 
some reason it is necessary to diminish the generated events and still have an umbiased 
sample for the whole process, this must be done choosing at random which events are 
to be taken. 

3 Program Structure 

In the following two sections we will explain how to use the input parameters to exploit 
the various possibilities of WPHACT, and the meaning of the various subroutines and 
functions of the code. At the end of the paper we report and briefly comment some 
significant examples of test runs. 

3.1 Input 

The way the input parameters must be given in the command file is easily under- 
stood just looking at the program lines containing the READ statements. We therefore 
reproduce them and explain their meaning in the following. 

Every parameter whose initial is i or n is integer*4. All others are real*8. When 
a variable has a yes/no option the value 1 corresponds to YES, to NO. All energies 
must be expressed in GeV. 

READ* , e_cni 

e_cni is the centre of mass energy. 

READ*,iproc 
READ*,ich 

The value for these two parameters must be read directly from Table 2 and 3. iproc 
selects the group of processes with the same cross section at parton level with symmetric 
cuts, while ich selects among them the specific final state. The choice of ich is relevant 
only if hadronization is performed via Jetset or if in CC processes cuts are considered 
for which a final configuration and its charge conjugate behave differently. In all other 
cases ich=l gives the same results as any other value indicated in the tables. 

IF ( iproc. GE. 33) THEN 
READ* , rmb 
READ*,icch 
IF(icch.NE.2)THEN 
READ*,isusy 
IF(isusy.EQ.O)THEN 
READ* , rmh 
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ELSE 
READ* , rma 
READ*,tgb 
ENDIF 
ENDIF 
ENDIF 

The processes for which iproc > 33 are those with massive 6's in the final state, rmb 
fixes the b mass for the phase space and the matrix elements. The b mass in the Higgs 
coupling may be different from rmb and it has to be set in the DATA. In these processes 
one can have diagrams with SM or Susy MSSM neutral {A or h) Higgs. One has the 
possibility of choosing to compute only these diagrams (Higgs signal: icch=l), only 
those without Higgs (to which we refer as Higgs background: icch=2) or the whole set 
of diagrams (Higgs+Background+interference: icch=3). 

isusy=0 corresponds to SM Higgs, isusy=l to Susy Higgs. In the first case one has 
to specify the Higgs mass (rmh), in the second the mass of the pseudoscalar Higgs A 
(rma), and the ratio between the two vacuum expectation values tg{P) (tgb). 

IF (iproc. GE. 6. AND. iproc. LE. 8) THEN 

READ* , iccnc 
ENDIF 

The processes for which 6 < iproc < 8 are those which have both CC and NC di- 
agrams (Mixed). For these one may choose to compute only the CC contribution 
(iccnc=l), the NC one (iccnc=2) or the total process including CC+NC+interference 
(iccnc=3). 

READ*,ips_cc 
READ* , ips_nc 

ips_cc and ips_nc allow to choose among the various phase space mappings for the 
integration. ips_cc refers to the phase space of CC or Higgs signal contributions. 
ips_nc to NC contributions. Of course when only CC processes or only Higgs signal 
are considered ips_nc is irrelevant. The same happens to ips_cc when only neutral 
processes without Higgs are computed. Both parameters can assume 3 values: 1 for 
double resonant mapping, 2 for single resonant, 3 for non resonant. When 2 is chosen 
the invariant mass over which a transformation is performed to take care of the resonant 
peak is that formed by particles 3 and 4. 

READ*,icc3 

yes/no CC3 contribution only. When icc3 is set to 1 only the three double resonant 
diagrams (CC3) corresponding to WW production and decay are computed. If icc3=0 
all CC contributions are computed. 

READ*,isr 

yes/no ISR: initial state radiation (ISR) is included when isr=l, not computed if isr=0. 
READ* , ipr 



12 



yes/no running widths: this flag selects among running or constant Z,W, Higgs widths 
in s-channel propagators: 

ipr=0 Z, W, Higgs boson constant width 

ipr=l Z, W, Higgs boson s-dependent width 

READ* , iswgcomp 

yes/no sin^Oy/ and g computed. If this flag is set to 1, it is used the renormahzation 
scheme in which siri?9w and g are computed from Z mass, W mass, Gf. If it is set to 
0, the values for sin^Ow and are taken from the DATA. The relation between g and 
ctem is always g'^ — Airaem/ sin^Ow 

READ* , igwcomp , igzcomp , ighcomp 

yes/no W, Z, H width computed. When one of these flags is =1, the corresponding 
W, Z or Higgs width is computed by standard formulas. If it is =0, the value for the 
corresponding width is the one given in the DATA. 

READ*,icoul 

Coulomb corrections may (icoul=l) or not (icoul=0) be computed. 
READ*,istrcor 

yes/no 'naive' QCD corrections to the cross sections. It has to be noticed that, if 
istrcor=0 and quarks are present in the flnal state, the eventual width computations 
(performed when igwcomp, igzcomp or ighcomp is =1) do not include QCD correc- 
tions. 

READ* , iqu 

yes/no QCD diagrams for 4-quarks NC. The order at which one computes four fermion 
diagrams is a^. There are some diagrams in NC processes O^aag) whose contribution 
can be important and can be included choosing iqu=l. 

READ*, i cut 
IF(icut.EQ.l)THEN 

READ* , e_min 

READ* , e_mcLX 

READ* , rm_min 

READ* , rm_max 

READ*,pt_min 

READ*,pt_max 

READ* , icos 

READ* , thbeam_min 

READ* , thbeam_max 

READ* , thsep_min 

READ* , thsep_max 
ENDIF 

Cuts may (icut=l) or not (icut=0) be implemented. If icut=l all default cuts of the 
above list must be specifled. 
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e_min and e_max correspond to the 4 lower and upper energies for particle 3, 4, 5, 6 
respectively. 

rmjuin and rmjnax are the 6 invariant mass lower and upper limits respectively. They 
must be given in the following order: m(34), m(35), m(36), m(45), m(46), m(56). 
pt_niin and pt_max are the 4 lower and upper values of the transverse momenta. The 
order is as before 3, 4, 5, 6. 

icos = 1 implies that the following angular cuts must be expressed in terms of the 
cosines of the angles. With icos = one must instead specify the angles in degrees. 
thbeam_niin and thbeam_max are the 4 lower and the 4 upper limits for the angle that 
particles 3,4,5,6 produce with the beam (e"*"). 

thsep_min and thsep_max are the 6 lower and the 6 upper limits for particle-particle 
angles: the order is again (3 4), (3 5), (3 6), (4 5), (4 6), (5 6). 

READ* , ianc 
IF(ianc.EQ.l)THEN 

READ*, delZjXf ,xz,yf ,yz,zz 
END IF 

Anomalous Couplings contributions may (ianc=l) or not (ianc=0) be computed. 
If one wants to compute them, he must also specify the values of the parameters 



5z, x^, xz, 2/7, yz, zz defined in ref. pT 



READ*,idistr 
IF(idistr.EQ.l)THEN 

READ*,ndistr 

DO i=l,ndistr 
READ*,nsubint(i) 

READ*, (distr_estrinf (i, j) , j=l,nsubint(i)+l) 
READ*, (nbin_number(i, j) , j=l,nsubint(i)) 
ENDDO 
END IF 

Distributions at parton level can be easily implemented. If one wants to use this 
possibility the flag idistr must be =1, =0 otherwise. In the program there is the line: 
* include 'abdis.dis' 

which must be uncommented before compiling if one choose idistr=l. In this case 
the file abdis.dis must be written for implementing distributions, as in the following 
example: 

string(l)='Distribution: bb~ invariant mass' 
distr.var ( 1) =sqrt ( (p3 (0) +p4 (0) ) **2- (p3 ( 1 ) +p4 ( 1) ) **2- 
& (p3(2)+p4(2))**2-(p3(3)+p4(3))**2) 

string(2)='Distribution: Charged lepton energy ' 
distr_var(2)=p5(0) 

In it, the title of the ith distribution is given in the character*60 string(i) and the 
ith quantity to be distributed in bins, computed from the particle momenta, is assigned 
to distr_var(i). The resulting cross sections corresponding to every single bin will be 
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stored in the file abdis .dat. Each fine will contain 3 numbers: the value of the central 
point of the bin, the distribution for the bin (cross section divided by the width of the 
bin, in order to reproduce da/dx for a distribution of the variable x) and the estimated 
statistical error. 

ndistr is the number of distributions defined in abdis. dis. 
For each distribution i, one must specify: 

nsubint(i), the number of sub- intervals with different binning in the ith distribu- 
tion (=1 when all bins are of the same length). The subintervals must be contiguous. 

(distr_estrinf (i, j) , j = l,nsubint(i)+l), the lower limits of each subinterval 
(which coincide with the upper limit of the previous one as they must be contiguous) 
and, as last entry, the upper limit of the last subinterval. In case all bins are of the 
same length, this corresponds only to the lower and upper limit of the interval for the 
distribution. 

(nbin_number(i, j) , j=l,nsubint(i)), the number of bins in each subinterval. 

READ*,iflat 
IFCiflat.EQ. DTHEN 

READ* , scalemax 

READ* , istorvegas 

READ* , irepeat 

IF ( irepeat . eq . 2) THEN 
READ*, nflevts 

END IF 

READ* , istormom 
READ*, i jet set 

IF (iproc . GE . 6 . AND . iproc . LE . 8) THEN 

READ*,interf 
ENDIF 
ENDIF 

One may choose (if lat=l) or not (if lat =0) to generate unweighted events. In the 
first case, the number of iterations (itmx) which must be specified in the following 
must be 2. The integration routine will perform the requested number of iterations for 
thermalization (see iterm, ncall_term, itmx_term below) and then the two iterations 
in which the integral is evaluated. In the first iteration the maximum for the hit-or-miss 
procedure will be determined and used in the second iteration where the unweighted 
generation will take place. After the run the .log file will report as usual the result of 
the integration and its error. It will also report the maximum used for the hit-or-miss 
procedure, the maximum found in the second iteration, and the number of events which 
were greater than the maximum used. 

There is also the possibihty to repeat the generation just starting directly from the 
second iteration. This might be useful if too many events exceeded the maximum 
chosen, or to generate a predetermined number of events. 

scalemajx is the coefficient by which the maximum of the first iteration can be multi- 
plied, in order to vary the efficiency of the hit-or-miss procedure or in order to avoid 
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values exceeding the maximum. 

VEGAS data are (if istorvegas=l) or not (if istorvegas=0) stored after the first 
iteration in ABVEGAS.DAT (in ABVEGAS_CC.DAT and ABVEGAS_NC.DAT for Mixed and 
Higgs+background processes). Stored VEGAS data are necessary if one wants to re- 
run the program to generate again unweighted events. When the program is rerun 
using VEGAS data stored, the maximum of the second iteration will be automatically 
used as the new maximum. 

irepeat has to be set to for the first run. It has to be set to 1 if one wants to rerun 
exactly with the same input starting from the second iteration. In this case the same 
weighted points will be reproduced. irepeat=2 has to be chosen if one wants to rerun 
with the same input and grid as before, but letting the program run until a requested 
number of events nf Itevts is reached. For both cases irepeat=l and 2 one might 
of course vary scalemax, ijetset and istormom with respect to the first run with 
irepeat=0. 

The momenta of the unweighted events are written in ABMOM.DAT (or ABMOM_SIGN.DAT 
and ABMOM_BACK.DAT) file if istormom=l, they are not written if istormom=0. 
ABMOM_SIGN.DAT is used for CC or Higgs events, ABMOM_BACK.DAT for NC events when 
mixed or Higgs+background processes are computed. 

Every unweighted event is passed to the standard COMMON HEPEVT. If ijetset=l the 
subroutine AB_LU4FRM is called for interfacing Jetset. If ijetset=0, HEPEVT is still filled 
but the interface to Jetset not called. 

One may choose to produce events for mixed processes in which the interference is added 
to CC (interf= 1), or to NC (interf=0). This implies that also the interference will 
be considered to have the same colour structure of CC or NC respectively (see section 
5.3 of ref. [|T§. 



READ*,acc 
READ*,iterm 
READ*,ncall_term 
READ* , itmx_term 
READ*,ncall 
READ* , itmx 

These parameters specify how the integration will be performed by VEGAS, 
acc is the integration accuracy. When this accuracy is reached after a certain integra- 
tion iteration, the remaining iterations are not performed. 

If iterni=l a certain number of integration iterations (=itmx_term) are used only for 
adapting the integration grid. Their result is not used for the final integral. Each ther- 
malizing iteration makes use of a maximum of ncall_term evaluations of the amplitude. 
If iterm=0 these iterations are skipped. One or two itmx_term with few ncall_teriii 
are often useful. 

itmx is the maximum number of iterations used to evaluate the integral. Each itera- 
tion makes use of a maximum of ncall evaluations of the amplitude. Normally it is 
better not to use more than about five iterations. If higher precision is requested it is 
convenient to increase ncall and not itmx 

As a final remark about the choice of these parameters, one must be aware of the fact 
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that final results with a much greater than the number of iterations are not to be 
trusted. When this happens, one has to increase ncall. 

In addition to the parameters in input, other parameters are fixed in the main 
program by the following DATA statement and may be eventually changed. 

DATA rmw/80 . 26d0/ , rmz/91 . 1884d0/, rmt/175.d0/, rmc/0 . 75d0/ , 
& rintau/1.78d0/, rmb_run/2 . 9d0/ , 
& gamw/2 . 08d0/ , gamz/2 . 4974d0/ , gainh/1.278d-03/, 

& gf/1.1663892199930875d-05/, alfainv/128.07d0/, 
& alfas_cc/0. 1255d0/, alf as_nc/0 . 1230d0/ , 
& s2w/0.231030912451068d0/, rms/lOOO.dO/ 

rmw, rmz, rmt, rmc, rmtau are respectively the VF, Z, top, c and tau masses. 
rmb_run is the quark b mass used for the Higgs coupling, 
gamw, gamz, gamh are the total W. Z and Higgs width. 

gf is the Fermi coupling constant; alfainv is l/ctem at the appropriate scale; 
alf as_cc and alf as_nc are as{Mw) and as{Mz)- 
s2w is the Weinberg sin^ {9w) and rms the SUSY scale. 

3.2 Program implementation 

In this section, WPHACT subroutines and functions are briefly described, 
double precision FUNCTION fxn(x,wgt) 

COMMON/ abpara/ rmxl , gamxl , gxl , rmx2 , gamx2 , gx2 , bet a , rlim, s_col , 
& xl_inin , x2_inin , xm , smin , emcoupl , estrinf 1 , estrinf 2 , estrined2 , 
& rmx3 , gamx3 , gx3 

COMMON/ abparb/ rmw2 , gamw , rmz2 , rmt2 , rcotw , pi , alf a_ine , qcdcoupl , 
& qcdcor_cc,qcdcor_nc 

COMMON/ abopzi/ isr , ipr , ips , iccnc , iproc , ich , ichc j , ips_cc , ips_nc , 
& icos , icch, isusy , icut , igwcomp, igzcomp, ighcomp, iswgcomp 

COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 
& f 3r , f 4r , f 5r , f 6r , z3r , z4r , z5r , z6r , wcl , delz , xf , xz , yf , yz , zz 

COMMON/ abf lag/ ice , icc3 , i3e , i4e , i3q, i5q, iqu , i34e , i56ve , ibbveve , iid, 
& imix , icoul , istrcor , idownl , idownr , ieinc 

COMMON/abhigg/rmb , riiib2 , rmh , rinh2 , gamh , rhzz , rhww , rhbb 

COMMON/absusy/rma, rnia2 , rzha 

C0MM0N/abcuts/e_min(3:6) ,e_max(3:6) ,thbeain_min(3 : 6) , 
& thbeam_max(3:6) ,thsep_min(6) ,thsep_max(6) ,rm_min(6) ,rm_max(6) , 
& beamcut_min(3 : 6) ,beamcut_mcLx(3 : 6) , sepcut_inin(6) , sepcut_max(6) , 
& rm_iiiin2 (6) ,rm_max2(6) ,pt_min(3:6) ,pt_inax(3:6) ,e_cm 

COMMON/abstat/ncall_eff 

COMMON/abdist/distr_estrinf (ndismaix.nestnnax) ,bin_width(ndismax, 
& nintmax) ,distr_local(ndismax,nbinmax,iiitmax) ,distr_loc_mix 
& (ndismax,nbinmax,nitiiiajc) ,dev_local(ndisinax,nbinmax,nitmEix) , 
& dev_loc_mix(ndismaj!;,nbinmax,nitmax) ,tail_local(ndismax,nitmax) , 
& tail_loc_mix(ndismcLx,nitmax) 
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COMMON/abidis/ncallbin(ndismeix,nbiiimax,nitmcLx) , 

& iicallbin_mix(ndisinax,nbininax,nitmax) ,nbin_nuinber(ndismax, 
& nintmax) ,nbin_suin(ndismax,nintinEix) ,nsubint(ndismax) .ndistr, 
& idistr,itl,it2,init 

COMMON/abcdis/stringCndismax) 

COMMON/ abf lat/rmaxf xn.rmaxf xn_lit ,rmaxf xn_2it ,rm£Lxf xn_cc_lit , 
& rmaxf xn_cc_2it .rmaxf xn_nc_lit ,rmaxf xn_nc_2it , scalemax 

COMMON/abif la/itmx ,noverinax, if lat , iseed, istorvegas , istormom, item, 
& ijetset , interf 

COMMON/abf la2/irepeat .nevent ,nf levts 

COMMON /abrann/ idum 

COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP) ,IDHEP(NMXHEP) , 
COMMON/HEPEVT/NEVHEP , NHEP , ISTHEP (NMXHEP) , IDHEP (NMXHEP) , 
& JMOHEP (2, NMXHEP ) ,JDAHEP( 2, NMXHEP) ,PHEP( 5, NMXHEP) ,VHEP( 4, NMXHEP) 
EXTERNAL ee_4f ,ee_bbvv,ee_bbmiimu,ee_bbbb,ee_bbee 

The function f xn is called by VEGAS and evaluates the phase space and ISR. After mo- 
menta reconstruction, it eventually implements the cuts and then calls the appropriate 
functions to compute the matrix elements. It also performs all calculations relative 
to distributions, it fills the common HEPEVT, and if necessary it calls the interface to 
Jetset, AB_LU4FRM. 

double precision FUNCTION ee_4f (pl,p2,p3,p4,p5,p6) 

COMMON/ abparb/ rmw2 , gamw , rmz2 , rmt2 , rcotw , pi , alf a_me , qcdcoupl , 

& qcdcor_cc,qcdcor_nc 
COMMON/ abpar c/ czipr , ccz , cwipr , ccw , chipr , cch , caipr , cca 
COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 

& f 3r , f 4r , f 5r , f 6r , z3r , z4r , z5r , z6r , wcl , delz , xf , xz , yf , yz , zz 

COMMON/abf lag/ice , icc3 , i3e , i4e , i3q, i5q, iqu , i34e , i56ve , ibbveve , iid, 

& imix , icoul , istrcor , idownl , idownr , ianc 

ee_4f computes the amplitude for all the processes with massless four fcrmions in the 
final state. It also evaluates Coulomb, QCD corrections and Anomalous Couplings. It 
is called by f xn. 

double precision FUNCTION ee_bbvv(pl,p2,p3,p4,p5,p6) 
COMMON/abhigg/ rmb , rmb2 , rmh , rmh2 , gamh , rhzz , rhww , rhbb 
C0MM0N/abparb/rmw2 , gamw , rinz2 , rmt2 , rcotw , pi , alf a_me , qcdcoupl , 

& qcdcor_cc,qcdcor_nc 
COMMON/ abpar c/ czipr , ccz , cwipr , ccw , chipr , cch , caipr , cca 
COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 

& f 3r,f4r,f5r ,f6r,z3r,z4r ,z5r,z6r , wcl , delz, xf ,xz,yf ,yz ,zz 
COMMON/ abcopl/ zvl , zvr , f qdl , f qdr , zqdl , zqdr , f qui , f qur , zqul , zqur 
COMMON/abf lag/ice , icc3 , i3e , i4e , i3q, iSq, iqu , i34e , i56ve , ibbveve , iid, 

& imix , icoul , istrcor , idownl , idownr , iainc 

ee_bbvv is called by f xn and gives the massive amplitude relative to the processes: 
e+e" — >• bbuiui l={e,ii,T). 

double precision FUNCTION ee_bbmumu(pl,p2,p3,p4,p5,p6) 
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COMMON/abhigg/rmb , rmb2 , rmh , rinh2 , gamh , rhzz , rhww , rhbb 
CDMMDN/abparb/rmw2 , gamw , rinz2 , rint2 , r cotw , pi , alf a_me , qcdcoupl , 
& qcdcor_cc,qcdcor_nc 
COMMON/abparc/ czipr , ccz , cwipr , ccw , chipr , cch , caipr , cca 

COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 
& f 3r ,f 4r ,f 5r ,f 6r ,z3r ,z4r ,z5r ,z6r ,wcl ,delz,xf ,xz,yf ,yz ,zz 
COMMON/ abcopl/zvl , zvr , f qdl , f qdr , zqdl , zqdr , f qui , f qur , zqul , zqur 
COMMON/abf lag/ice , icc3 , i3e , i4e , i3q, i5q, iqu, i34e , i56ve , ibbveve , iid, 
& imix , icoul , istrcor , idownl , idownr , ianc 

ee_bbmumu is called by fxn and computes the massive amplitude relative to the pro- 
cesses: 

e^e~ —>■ bbqq q={u,d,c,s) 
e+e-^bbl-l+ l^{f^,T) 

double precision function ee_bbee(ql,q2,q3,q4,q5,q6) 

COMMON/abhigg/rmb , rmb2 , rmh , rmh2 , gamh , rhzz , rhww , rhbb 

COMMON/ abparb/ rmw2 , gamw , rmz2 , rmt2 , r cotw , pi , alf a_me , qcdcoupl , 

& qcdcor_cc,qcdcor_nc 
COMMON/abparc/czipr , ccz , cwipr , ccw , chipr , cch , caipr , cca 
COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 

& f 3r ,f 4r ,f 5r ,f 6r ,z3r ,z4r ,z5r , z6r ,wcl ,delz ,xf ,xz ,yf ,yz , zz 
COMMON/ abcopl/ zvl , zvr , f qdl , f qdr , zqdl , zqdr , f qui , f qur , zqul , zqur 
COMMON/abf lag/ ice , icc3, i3e , i4e , i3q, i5q, iqu, i34e , i56ve , ibbveve , iid, 

& imix , icoul , istrcor , idownl , idownr , ianc 

ee_bbee is called by fxn and gives the massive amplitude relative to the process : 
e'^e~ — > bbe~e'^ 

double precision FUNCTION ee_bbbb(pl ,p2,p3,p4,p5,p6) 

COMMON/ abopzi/ isr , ipr , ips , iccnc , iproc , ich , ichc j , ips_cc , ips_nc , 
& icos,icch, isusy , icut , igwcomp , igzcomp , ighcomp , iswgcomp 

C0MM0N/abparb/rmw2 , gamw , rmz2 , rmt 2 , r cotw , pi , alf a_me , qcdcoupl , 
& qcdcor_cc,qcdcor_nc 

COMMON/ abparc/ czipr , ccz , cwipr , ccw , chipr , cch , caipr , cca 

COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 
& f 3r , f 4r , f 5r , f 6r , z3r , z4r , z5r , z6r , wcl , delz , xf , xz , yf , yz , zz 

COMMON/ abcopl/ zvl , zvr , f qdl , f qdr , zqdl , zqdr , f qui , f qur , zqul , zqur 

COMMON/abf lag/ice , icc3 , i3e , i4e , i3q, iSq, iqu , i34e , i56ve , ibbveve , iid, 
& imix , icoul , istrcor , idownl , idownr , ianc 

COMMON/abhigg/rmb , rmb2 , rmh , rmh2 , gamh , rhzz , rhww , rhbb 

COMMON/absusy/rma, rma2 , rzha 

ee_bbbb computes the massive amplitude relative to the process e'^e~ — > bbbb . It 
contains both SM (hZ) and MSSM {hZ + hA) Higgs diagrams. It is called by fxn. 

Subr out ine AB_LU4FRM ( i char , TRAD , ITAU , lERR) 
C0MM0N/LUJETS/N,K(4000,5) ,P(4000,5) ,V(4000,5) 
C0MM0N/LUDAT1/MSTU(200) ,PARU(200) ,MSTJ(200) ,PARJ(200) 
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This routine is an interface to Jetset that we have produced modifying the subroutine 
LU4FRM by T. Sjostrand to adapt it to WPHACT. The comment hues of the routine 
AB_LU4FRM may be useful 

SUBROUTINE vegas (r egion , ndim , f xn , init , ncall , itmx , nprn , tgral , sd , 
* chi2a, acc ,xi , it ,ndo , si , swgt , schi) 
EXTERNAL fxn 

CDMMDN/abresl/resKlO) .standdevKlO) 
COMMON /abrann/ idum 
COMMON/abchia/calls 
COMMON/abstat/iicall_ef f 
COMMON/abf la2/irepeat ,nevent ,nf levts 

This subroutine by P. Lepage performs the multidimensional integrations. It has 
been modified in order to obtain unweighted event generation, distributions and sepa- 
rated integration of CC ( or Higgs) and NC contributions both for MIX and massive b 
processes. It is called by the main program and it calls the function fxn. 

FUNCTION ran2(idum) 
COMMON/absalv/iv , iy , idum2 

This function is a random number generator. Both VEGAS and WPHACT make use of it. 
SUBROUTINE rebiii(rc,nd,r,xin,xi) 

This routine is used by VEGAS. 

double precision FUNCTION gainmln(xx) 
It computes the logarithm of the T function entering in the electron structure functions. 



SUBROUTINE initialize 

COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 
& f 3r ,f 4r ,f 5r ,f 6r ,z3r ,z4r ,z5r ,z6r ,wcl ,delz,xf ,xz,yf ,yz ,zz 

COMMON/ abcopl/ zvl , zvr , f qdl , f qdr , zqdl , zqdr , f qui , f qur , zqul , zqur 

COMMON/ abf lag/ ice , icc3 , i3e , i4e , i3q, i5q, iqu, i34e , i56ve , ibbveve , iid, 
& imix , icoul , istrcor , idownl , idownr , ianc 

COMMON/ abf lat/rmaxf xn.rmaxf xn_lit ,rmaxf xn_2it ,rmaxf xn_cc_lit , 
& rmaxf xn_cc_2it ,rmcLxf xn_nc_lit ,rmaxf xn_nc_2it , scalemax 

COMMON/ abif la/ itmx,novermax, if lat , iseed, istorvegas , istormom, iterm, 
& ijetset , interf 

COMMON/abf la2/irepeat .nevent ,nf levts 

COMMON/ abopzi/ isr , ipr , ips , iccnc , iproc , ich, ichc j , ips_cc , ips_nc , 
& icos , icch , isusy , icut , igwcomp , igzcomp , ighcomp , iswgcomp 

COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP) ,IDHEP(NMXHEP) , 
& JM0HEP(2,NMXHEP) , JDAHEP(2,NMXHEP) ,PHEP (5 , NMXHEP) ,VHEP(4,NMXHEP) 
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For a process chosen by means of input parameters, this subroutine, called by the main 
program, computes all chiral coupling constants and it initializes all variables to be 
passed to other parts of the program and to the common HEPEVT. 

SUBROUTINE pr inter (rmh , gamh , rma , gama , tgb , rmb) 

COMMON/ abopzi/ isr , ipr , ips , iccnc , iproc , ich , ichc j , ips_cc , ips_nc , 
& icos , icch, isusy , icut , igwcomp, igzcomp, ighcomp, iswgcomp 

COMMON/ abparb/ rmw2 , gamw , rmz2 , rmt2 , r cotw , pi , alf a_me , qcdcoupl , 
& qcdcor_cc,qcdcor_nc 

COMMON/abcoup/f er , f el , zer , zel , f 31 , f 41 , f 51 , f 61 , z31 , z41 , z51 , z61 , 
& f 3r ,f 4r ,f 5r ,f 6r ,z3r ,z4r ,z5r ,z6r ,wcl ,delz,xf ,xz,yf ,yz ,zz 

COMMON/ abf lag/ice , icc3 , i3e , i4e , i3q, i5q, iqu, i34e , i56ve , ibbveve , iid, 
& imix , icoul , istrcor , idownl , idownr , ianc 

C0MM0N/abcuts/e_min(3:6) ,e_max(3:6) ,thbeain_min(3 : 6) , 
& thbeajn_max(3 : 6) ,thsep_min(6) ,thsep_max(6) ,rm_min(6) ,rm_max(6) , 
& beajncut_min(3 : 6) ,beamcut_max(3:6) , sepcut_min(6) , sepcut_max(6) , 
& rm_miii2 (6) ,rm_max2(6) ,pt_min(3 : 6) ,pt_max(3:6) ,e_cm 

COMMON/ abinpu/ rmw , rmz , rmb_run , rmc , rnic_run , rmt au , gamz , gf , s2w , 
& alfaiiiv,alfas_cc,alfas_nc 

COMMON/ abf lat/rmaxf xn.rmaxf xn_lit ,rmaxf xn_2it ,rma:xf xn_cc_lit , 
& rmaxf xn_cc_2it ,rmaxf xn_iic_lit ,rmaxf xn_nc_2it , scalemax 

COMMON/ abif la/ itmx,novermax, if lat , iseed, istorvegas , istormom, iterm, 
& ijetset , interf 

COMMON/abf la2/irepeat .nevent ,nf levts 

This subroutine is devoted to write in the output file all the essential informations 
about the selected process, input, options and variables defined in the DATA. It is called 
by the main program. 



4 Conclusions 

We have described version 1.0 of WPHACT . The program computes all processes with 
a 4 fermion final state. All its main features and options are listed in the Program 



Summary. The way of using it is mainly explained in the input section ^]T|, and some 
useful examples can be found in the test runs. 

WPHACT is a MC integrator and event generator portable and self contained. For 
the amplitudes computation it makes use of a rather new helicity method which allows 
fast evaluations. This we consider to be an essential ingredient, together with careful 
mappings of the integration variables, in order to combine the usual advantages of MC's 
with high precision and reliability. 
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Test runs 



We report here two significant test runs for CC and NC+Higgs processes witli different 
options. 

Run 1: 

Tliis run computes tfie cross section and two distributions for tlie full CC20 process 
e+e~ — > e'lygud. A double resonant phase space on the W masses has been chosen to 
perform the integration over the 9 variables. ISR, Coulomb and QCD corrections have 
been included as well as Anomalous Couplings. The following cuts arc implemented: 
Ee-,u4 > 20 GeV, M„,- > 10 GeV, \cos{e,,am)\e- ,u,d < 0.9, \cos{e^i)\ < 0.9. 
Distributions of the Mw+ invariant mass and the outgoing electron energy are calcu- 
lated. The corresponding plots are reported in figs. 1,2. They contain 51 and 150 bins 
respectively. The relative statistical errors are a few permill for every bin. They are 
not visible in the figures. 

The relative error obtained is 0.019%. The total CPU time on AlphaStation 600 5/333 
is 24 min and 36 sec, which corresponds to 0.6E — 04 sec. per call. 

Input file: 



$r wphact 
190. dO 
4 








1 

20. dO O.dO 20. dO 20. dO 

300. dO 300. dO 300. dO 300. dO 
O.dO O.dO O.dO O.dO O.dO 10. dO 



centre of mass energy 

selects the kind of process 

selects the channel 

CC( or Higgs signal) phase space 

NC phase space 

yes /no CCS only 

yes/no ISR 

yes/no running widths 

yes/no s2w and g computed 

yes/no W, Z, Higgs width computed 

yes/no Coulomb corrections 

yes /no QCD corrections 

yes/no QCD diagrams 

yes/no cuts 

4 energy lower cuts 

4 energy upper cuts 

6 invariant mass lower limits 



300. dO 300. dO 300. dO 300. dO 300. dO 300. dO !6 invariant mass upper limits 



O.dO O.dO O.dO O.dO 
300. dO 300. dO 300. dO 300. dO 
1 

0. 9d0 l.dO 0.9d0 0.9d0 
-0.9d0 -l.dO -0.9d0 -0.9d0 

1. dO l.dO l.dO l.dO l.dO 0.9d0 
-l.dO -l.dO -l.dO -l.dO -l.dO -0.9d0 
1 

-0.5d0 O.dO O.dO O.dO O.dO O.dO 
1 
2 
3 



4 transverse momenta lower cuts 
4 transverse momenta upper cuts 
angular cuts in deg (0) or cos (1) 
4 particle-beam angle lower cuts 
4 particle-beaim angle upper cuts 
6 particle-particle angle lower cuts 
!6 part-part angle upper cuts 
yes/no anomalous couplings 
cinomalous couplings parameters 
yes/no distributions 
number of distributions 

sub-intervals number with different binning 
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75.71d0 78.21d0 82.31d0 84.81d0 

5 41 5 

1 

15. dO 90. dO 

150 



O.OOOldO 

1 

2000000 
2 

10000000 
2 

$exit 



limits of each sub-interval 
number of bins in each sub-interval 
sub-intervals number with different b' 
limits of each sub-interval 
number of bins in each sub-interval 
yes/no flat generation 
integration accuracy 
yes/no thermalization 
thermalization calls per iteration 
thermalization iterations 
integration calls per iteration 
integration iterations 



File abdis.dis : 

string(l)= 'Distribution: Mw+' 

distr_var ( 1 ) =sqrt ( (p5 (0) +p6 (0) ) **2- (p5 ( 1) +p6 ( 1) ) **2- 
& (p5 (2) +p6 (2) ) **2- (p5 (3) +p6 (3) ) **2) 

string(2)='Distribution: Lepton energy' 
distr_var(2)=p3(0) 

Output file: 



CC20 ) e-(p3) ve~(p4) u(p5) d~(p6) 
WW signal + background 

INPUT 



cm energy 







1900000D+03 


GeV 


DATA 










Z mass 







91188400+02 


GeV 


W mass 







8026000D+02 


GeV 


Z width 







2497400D+01 


GeV 


Gf 







1166389D-04 


GeV-2 


s2w 







2310309D+00 




1/alf a_em 







1280700D+03 




alf as_cc 







1255000D+00 




DERIVED QUANTITIES 










W width 







2088612D+01 


GeV 



OPTIONS 

both Z and W boson s-dependent width 
Born + QED 

Coulomb corrections included 
Naive QCD corrections included 
Double resonant "CC" phase space 

Anomalous coupling parameters: 
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delz = 


-0 


5 


xf = 








xz = 








yf = 








yz = 








zz = 









Cuts : 



ENERGY_MIN(3,4,5,6) =( 


20.00, 


0.00, 20.00, 20 


00 ) 


GeV 






ENERGY_MAX(3,4,5,6) =( 


300.00, 


300.00, 300.00, 


300. 


00 ) GeV 






MASS_MIN(34,35,36,45,46,56) =( 


0.00, 


0.00, 0.00, 


00, 


0.00, 10 


00 


) GeV 


MASS_MAX(34,35,36,45,46,56) =( 


300.00, 


300.00, 300.00, 


300. 


00, 300.00, 300.) 


PT_MIN(3,4,5,6) =( 


0.00, 


0.00, 0.00, 


00 ) 


GeV 






PT_MAX(3,4,5,6) =( 


300.00, 


300.00, 300.00, 


300. 


00 ) GeV 






C0SBEAM_MAX(3,4,5,6) =( 


0.90, 


1.00, 0.90, 


90 ) 








C0SBEAM_MIN(3,4,5,6) =( 


-0.90, 


-1.00, -0.90, -0 


90 ) 








C0SSEP_MAX(34,35,36,45,46,56) =( 


1.00, 


1.00, 1.00, 1 


00, 


1.00, 


90 


) 


C0SSEP_MIN(34,35,36,45,46,56) =( 


-1.00, 


-1.00, -1.00, -1 


00, 


-1.00, -0 


90 


) 



Thermalization 

input parameters for vegas: 



ndiin= 9 ncall= 1835008. 
it= 1 itmx= 2 



iteration no . 1 : 
iteration no. 1: integral 
all iterations: integral 

iteration no. 2: 
iteration no. 2: integral 
all iterations: integral 

CC process 

input parameters for vegas: 



effective ncall 
0.4101226 
0.4101226 



1052601 
+/- 0.50E-03 

+/- 0.503E-03 chi**2/it'n = 0.12E-07 



effective ncall= 1273154 
0.4100550 +/- 0.29E-03 

0.4100716 +/- 0.249E-03 chi**2/it'n = 0.13E-01 



ndim= 9 ncall= 9765625. 
it= 1 itmx= 2 



iteration no . 1 : 

iteration no. 1: integral 

all iterations: integral 

iteration no. 2: 

iteration no. 2: integral 

all iterations: integral 



effective ncall 
0.4101164 
0.4101164 



6914065 
+/- O.llE-03 

+/- O.llOE-03 chi**2/it'n = O.OOE+00 



effective ncall= 6934162 
0.4100655 +/- O.llE-03 

0.4100908 +/- 0.774E-04 chi**2/it'n =0.11 



Sigma = 0.4100908D+00 +/- 0.774D-04 (pb) 



Run 2: 

This output refers to NC25 process e+e" —>■ bbii~ ii"^ with massive 6's in the final state. 
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The cross section receives contributions both from Higgs and Background NC diagrams 
integrated by means of two phase spaces accounting for ZH and ZZ peak structure 
respectively. ISR and QCD corrections are included. The cuts are: M^j > 50 GeV and 
Mz — 25 GeV < M^^ < Mz + 25 GeV. This output shows in particular two possible 
examples of unweighted events generation: the first run produces unweighted events 
(which are stored in ABMOM_ SIGN (BACK) .DAT) with an efficiency of a few percent. 
Making use of the Vegas data saved in VEGAS_ CC(NC) .DAT during the first run, a 
second run allows to generate a prefixed number of 50000 unweighted events in a short 
CPU time. 

The relative error obtained is 0.0028%. The total CPU time for the both runs on 
AlphaStation 600 5/333 is 4 h 4 min and 5 sec. 

First input file: 
$r wphact 

190. dO ! centre of mass energy 

36 ! selects the kind of process 

1 ! selects the channel 

2.9d0 ! quark "b" mass 

3 ! icch (=1 Higgs, =2 Backg. , =3 Higgs+Backg. ] 

! yes/no SUSY 
80. dO ! Higgs mass 

1 ! CC( or Higgs) phase space 
1 ! NC phase space 

! yes /no CC3 only 

1 ! yes/no ISR 

1 ! yes/no running widths 

1 ! yes/no s2w and g computed 

10 1 ! yes/no W, Z, Higgs width computed 

! yes/no Coulomb corrections 

1 ! yes/no QCD corrections 

! yes/no QCD diagrams 

1 ! yes/no cuts 

O.dO O.dO O.dO O.dO ! 4 energy lower cuts 

300. dO 300. dO 300. dO 300. dO ! 4 energy upper cuts 

50. dO O.dO O.dO O.dO O.dO 66.1888d0 ! 6 invariant mass lower limits 

300. dO 300. dO 300. dO 300. dO 300. dO 116.1888d0 !6 invariant mass upper limits 



O.dO O.dO O.dO O.dO 

300. dO 300. dO 300. dO 300. dO 



O.dO O.dO O.dO O.dO 

180. dO 180. dO 180. dO 180. dO 

O.dO O.dO O.dO O.dO O.dO O.dO 



180. dO 180. dO 180. dO 180. dO 180. dO 180. dO !6 part -part angle upper cuts 



IdO 



4 transverse momenta lower cuts 
4 transverse momenta upper cuts 
angular cuts in deg (0) or cos (1) 
4 particle-beam angle lower cuts 
4 particle-beam angle upper cuts 
6 particle-particle angle lower cuts 



yes/no anomalous couplings 
yes/no distributions 
yes/no flat generation 
scale factor for the maximum 
yes/no data VEGAS stored 
yes/no second iteration repeated 
yes/no flat momenta stored 
yes/no Jetset program tied 



27 



O.OOOSdO ! integration accuracy 

1 ! yes/no thermalization 

1000000 ! thermalization calls per iteration 

2 ! thermalization iterations 
20000000 ! integration calls per iteration 
2 ! integration iterations 

$exit 



First output file: 



NC25 ) b(p3) b~(p4) mu-(p5) mu+(p6) 
Higgs signal + background 



INPUT 
cm energy 
Higgs mass 
b mass 



0.1900000D+03 GeV 
0.8000000D+02 GeV 
0.2900000D+01 GeV 



DATA 



Z mass 







9118840D+02 


GeV 


W mass 







8026000D+02 


GeV 


c mass 







7500000D+00 


GeV 


tau mass 







1780000D+01 


GeV 


Z width 







2497400D+01 


GeV 


Gf 







1166389D-04 


GeV-2 


alf as_nc 







1230000D+00 





DERIVED QUANTITIES 
W width 
Higgs width 
s2w 

1/alf a_em 



0.2090172D+01 GeV 
0. 19373740-02 GeV 
0.2253258D+00 
0.1312146D+03 



OPTIONS 

both Z and H boson s-dependent width 
Born + QED 

Naive QCD corrections included 
Double resoneint "Higgs signal" phase space 
Double resonant "Higgs background" phase space 
Cuts : 



ENERGY_MIN(3,4,5,6) 


= ( 


0.00, 


ENERGY_MAX(3,4,5,6) 


= ( 


300.00 


MASS_MIN (34 , 35 , 36 , 45 , 46 , 56) 


= ( 


50.00, 


MASS_MAX (34 , 35 , 36 , 45 , 46 , 56) 


= ( 


300.00 


PT_MIN(3,4,5,6) 


= ( 


0.00, 


PT_MAX(3,4,5,6) 


= ( 


300.00 


THBEAM_MIN(3,4,5,6) 


= ( 


0.00, 


THBEAM_MAX(3,4,5,6) 


= ( 


180.00, 


THSEP_MIN (34 , 35 , 36 , 45 , 46 , 56) 


= ( 


0.00, 



0.00, 0.00, 0.00 ) GeV 
, 300.00, 300.00, 300.00 ) GeV 

0.00, 0.00, 0.00, 0.00, 66.19 ) GeV 
, 300.00, 300.00, 300.00, 300.00, 116.) 

0.00, 0.00, 0.00 ) GeV 
, 300.00, 300.00, 300.00 ) GeV 

0.00, 0.00, 0.00 ) deg 
180.00,180.00,180.00 ) deg 

0.00, 0.00, 0.00, 0.00, 0.00 ) deg 
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THSEP_MAX(34,35,36,45,46,56) =( 180.00,180.00,180.00,180.00,180.00,180.00 ) deg 



Flat events generation 

VEGAS data stored in ABVEGAS_CC(NC) .DAT 
Maximum scale factor = O.llOD+01 
Flat events stored in ABMOM_SIGN.DAT 
Flat events stored in ABMOM_BACK.DAT 

Thermalization 



input parameters for vegas : ndim= 9 ncall= 

it= 1 itmx= 



786432. 



iteration no . 1 : 
iteration no. 1: integral 
all iterations: integral 



effective ncall= 
0.1661491E-01+/- 
0.1661491E-01+/- 



743476 
0.75E-05 

0.745E-05 chi**2/it'n = O.OOE+00 



input parameters for vegas: 



ndim= 9 ncall= 
it= 1 itmx= 



786432 . 



iteration no . 1 : 
iteration no. 1: integral 
all iterations: integral 



effective ncall= 
0.7878380E-02+/- 
0.7878380E-02+/- 



743476 
0.61E-05 

0.608E-05 chi**2/it'n = O.OOE+00 



iteration no . 2 : 
iteration no. 2: 
all iterations: 



integral 
integral 



effective ncall= 
0.7882914E-02+/- 
0.7881695E-02+/- 



776518 
0.37E-05 
0.315E-05 chi**2/it'n =0.40 



Higgs signal 

input parameters for vegas: 



ndim= 9 ncall= 
it= 1 itmx= 



19531250. 
1 



iteration no . 1 : 
iteration no. 1: integral 
all iterations: integral 



effective ncall= 
0.1661659E-01+/- 
0.1661659E-01+/- 



19292845 
0.78E-06 

0.784E-06 chi**2/it'n = O.OOE+00 



Higgs backg + Higgs sign-backg interference 



input parameters for vegas: 



ndim= 9 ncall= 
it= 1 itmx= 



19531250. 
1 



iteration no . 1 : 
iteration no. 1: integral 
all iterations: integral 

Higgs signal 

input parameters for vegas: 



effective ncall= 
0.7879024E-02+/- 
0.7879024E-02+/- 



ndim= 9 ncall= 
it= 2 itmx= 



19345150 
0.62E-06 

0.621E-06 chi**2/it'n = O.OOE+00 



19531250. 
2 



iteration no. 2: 
iteration no. 2: integral 



effective ncall= 
0.1661579E-01+/- 



19317393 
0.76E-06 
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all iterations: integral = . 1661618E-01+/- 0.545E-06 chi**2/it'n =0.54 



Higgs backg + Higgs sign-backg interference 



input parameters for vegas 



nd.im= 9 ncall= 
it= 2 itmx= 



19531250. 
2 



iteration no . 



2: 



effective ncall= 



19355281 



iteration no. 2: integral = . 7877976E-02+/- 0.60E-06 

all iterations: integral = 0.7878485E-02+/- 0.433E-06 chi**2/it'n = 1.5 



Sigma = 0.2449466D-01 +/- 0.696D-06 (pb) 
Informations about flat events generation: 



Maximum after first VEGAS iteration = 0.309D-07 
Maximum after second VEGAS iteration = 0.215D-07 
Flat events number = 720755 

number of function values over maximum = 



Second input file: 

We report only the input changed with respect to the first input. 

l.dO ! scale factor for the maximum 

! yes/no data VEGAS stored 

2 ! yes/no second iteration repeated 

50000 ! number of unweighted events generated 

1 ! yes/no flat momenta stored 

Second output file: 

We do not report here INPUT, DATA, DERIVED QUANTITIES, OPTIONS and cuts which are 
the same as in the first output. 



MC25 ) b(p3) b~(p4) mu-(p5) mu+(p6) 
Higgs signal + background 



Flat events generation 
Maximum scale factor = O.lOOD+01 
Flat events stored in ABMOM_SIGN.DAT 
Flat events stored in ABMOM_BACK.DAT 

Higgs signal 

Flat events number = 33917 
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Higgs backg + Higgs sign-backg interference 
Flat events number = 16083 



Sigma = 0.2449466D-01 +/- 0.696D-06 (pb) 
Informations about flat events generation: 



Maximum = 1 . 330188364059000E-008 
Flat events number = 50000 

number of function values over maximum = 
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Table Caption 

Table 1 DEC AlphaStation 600 5/333 CPU time, accuracy and effective calls (in mil- 
lions) for some representative four fermion processes with ISR. 

Table 2 Charged Current and Mixed Charged + Neutral Current four fermion pro- 
cesses, iproc and ich are the two flags by which the appropriate final state is 
singled out in WPHACT. 

Table 3 Neutral Current and Neutral Current + Higgs four fermion processes, iproc 
and ich are the two flags by which the appropriate final state is singled out in 
WPHACT. 
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Table 1 
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Table 2 
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Table 3 
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Figure 1: Distribution of the invariant mass corresponding to M^+, computed in test 
run 1. 
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Figure 2: Distribution of the energy of the outgoing electron computed in test run 1. 
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